home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dlngam.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  2.0 KB  |  56 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((xmax 0.0)
  12.       (dxrel 0.0)
  13.       (sq2pil 0.9189385332046728)
  14.       (sqpi2l 0.22579135264472744)
  15.       (pi_ 3.141592653589793)
  16.       (first nil))
  17.   (declare (type f2cl-lib:logical first)
  18.            (type double-float pi_ sqpi2l sq2pil dxrel xmax))
  19.   (setq first f2cl-lib:%true%)
  20.   (defun dlngam (x)
  21.     (declare (type double-float x))
  22.     (prog ((sinpiy 0.0) (y 0.0) (temp 0.0) (dlngam 0.0))
  23.       (declare (type double-float dlngam temp y sinpiy))
  24.       (cond
  25.        (first (setf temp (/ 1.0 (f2cl-lib:flog (f2cl-lib:d1mach 2))))
  26.               (setf xmax (* temp (f2cl-lib:d1mach 2)))
  27.               (setf dxrel (f2cl-lib:fsqrt (f2cl-lib:d1mach 4)))))
  28.       (setf first f2cl-lib:%false%)
  29.       (setf y (coerce (abs x) 'double-float))
  30.       (if (> y 10.0) (go label20))
  31.       (setf dlngam (coerce (f2cl-lib:flog (abs (dgamma x))) 'double-float))
  32.       (go end_label)
  33.      label20
  34.       (if (> y xmax)
  35.           (xermsg "SLATEC" "DLNGAM" "ABS(X) SO BIG DLNGAM OVERFLOWS" 2 2))
  36.       (if (> x 0.0)
  37.           (setf dlngam
  38.                   (+ (- (+ sq2pil (* (- x 0.5) (f2cl-lib:flog x))) x)
  39.                      (d9lgmc y))))
  40.       (if (> x 0.0) (go end_label))
  41.       (setf sinpiy (coerce (abs (sin (* pi_ y))) 'double-float))
  42.       (if (= sinpiy 0.0)
  43.           (xermsg "SLATEC" "DLNGAM" "X IS A NEGATIVE INTEGER" 3 2))
  44.       (if (< (abs (/ (- x (f2cl-lib:aint (- x 0.5))) x)) dxrel)
  45.           (xermsg "SLATEC" "DLNGAM"
  46.            "ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER" 1 1))
  47.       (setf dlngam
  48.               (- (+ sqpi2l (* (- x 0.5) (f2cl-lib:flog y)))
  49.                  x
  50.                  (f2cl-lib:flog sinpiy)
  51.                  (d9lgmc y)))
  52.       (go end_label)
  53.      end_label
  54.       (return (values dlngam nil)))))
  55.  
  56.